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Abstract 

We analyze the spectrum of the ultrahigh energy (above « 
10 9 GeV) cosmic rays. With a maximum likelihood analysis 
we show that the observed spectrum is consistent with the 
decay of extragalactic GUT scale particles. The predicted 
mass for these superheavy particles is mx = 10 fc GeV, where 
6 = 14.6±i;f. 

1 Introduction 

The interaction of protons with photons of the cosmic mi- 
crowave background radiation (CMBR) predicts a sharp drop 
in the cosmic ray flux above the Greisen-Zatsepin-Kuzmin 
(GZK) cutoff around 5 • 10 19 eV (Greisen, 1966; Zatsepin 
and Kuzmin, 1966). The available data shows no such drop. 
About 20 events above 10 20 eV were observed by experi- 
ments such as AGASA (Takeda et al., 1998), Fly's Eye (Bird 
et al., 1993), Haverah Park (Lawrence et al., 1991), Yakutsk 
(Efimov et al., 1991) and HiRes (Kieda et al., 1999). Fu- 
ture experiments, particularly Pierre Auger (Boratav, 1996; 
Guerard, 1999; Bertou et al. 2000), will have a much higher 
statistics. 

Usually it is assumed that at these energies the galactic 
and extragalactic (EG) magnetic fields do not affect the orbit 
of the cosmic rays, thus they should point back to their ori- 
gin within a few degrees. Though there are clustered events 
(Hayashida et al., 1996; Uchihori et al., 2000) the overall 
distribution is practically isotropic (Dubovski and Tinyakov, 
1998; Berezinsky and Mikhailov, 1999; Medina Tanco and 
Watson, 1999), which usually ought to be interpreted as a 
signature for EG origin. Since above the GZK energy the 
attenuation length of particles is a few tens of megaparsecs 
(Yoshida and Teshima, 1993 ; Aharonian and Cronin, 1994; 
Protheroe and Johnson, 1996; Bhattacharjee and Sigl, 2000; 
Achterberg et al., 1999; T. Stanev et al., 2000) if an ultrahigh 
energy cosmic ray (UHECR) is observed on Earth it must be 
produced in our vicinity. Sources of EG origin should result 
in a GZK cutoff, which is in disagreement with experiments. 
It is generally believed that there is no conventional astro- 
physical explanation for the observed UHECR spectrum. 

An interesting idea is that superheavy particles (SP) as 
dark matter could be the source of UHECRs. In (Kuzmin and 
Rubakov, 1998) EG SPs were studied. A crucial observation 



was made (Berezinsky et al., 1997) about the decay of SPs 
concentrated in the halo of our galaxy. They used the modi- 
fied leading logarithmic approximation (MLLA) (Azimov et 
al., 1985; Fong and Webber 1991) for ordinary QCD and for 
supersymmetric QCD (Berezinsky and KachelrieB, 1998). A 
good agreement of the EG spectrum with observations was 
noticed in (Berezinsky et al., 1998). Supersymmetric QCD 
is treated as the strong regime of the minimal supersymmet- 
ric standard model (MSSM). To describe the decay spectrum 
more accurately HERWIG Monte-Carlo (Marchesini et al., 
1992) was used in QCD (Birkel and Sarkar, 1998) and dis- 
cussed in supersymmetric QCD (Rubin 1999, Sarkar 2000), 
resulting in m x « 10 12 GeV and « 10 13 GeV for the SP 
mass in SM and in MSSM, respectively. 

SPs are very efficiently produced by the various mecha- 
nisms at post inflatory epochs (for a review see Berezinsky 
2000). Note, that any analysis of SP decay covers a much 
broader class of possible sources. Several non-conventional 
UHECR sources produce the same UHECR spectra as de- 
caying SPs. 

Here we study the scenario that the UHECRs are com- 
ing from decaying SPs and we determine the mass of this X 
particle mx by a detailed analysis of the observed UHECR 
spectrum. We discuss both possibilities that the UHECR pro- 
tons are produced in the halo of our galaxy and that they are 
of EG origin and their propagation is affected by CMBR. 
Here we do not investigate how can they be of halo or EG 
origin, we just analyze their effect on the observed spectrum 
instead. We assume that the SP decays into two quarks (other 
decay modes would increase mx in our conclusion). Af- 
ter hadronization these quarks yield protons. The result is 
characterized by the fragmentation function (FF) D(x,Q 2 ) 
which gives the number of produced protons with momen- 
tum fraction x at energy scale Q. For the proton's FF present 
accelerator energies can be used (Binnenwies et al., 1995; 
Kniehl et al., 2000). We evolve the FFs in ordinary (Gri- 
bov and Lipatov 1972; Lipatov 1975; Altarelli and Parisi, 
1977; Dokshitzer, 1977) and in supersymmetric (Jones and 
Llewellyn Smith 1983) QCD to the energies of the SPs. This 
result can be combined with the prediction of the MLLA 
technique , which gives the initial spectrum of UHECRs at 
the energy m x . Altogether we study four different models: 
halo-SM, halo-MSSM, EG-SM and EG-MSSM. 
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Table 1 : The fragmentation functions of the different partons 
using the parametrization D{x) = Nx a (l — x) 13 at different 
energy scales (second column). 

2 Decay and fragmentation of heavy 
particles 

The UHECRs are most likely to be dominated by protons 
(Dawson et al., 1998); in our analysis we use them exclu- 
sively. 

The FF of the proton can be determined from present ex- 
periments (Binnenwies et al., 1995; Kniehl et al., 2000). The 
FFs at Qq energy scale are Di(x, Qq), where i represents the 
different partons (quark, squark or gluon, gluino). The FFs 
can not be determined in perturbative QCD; however, their 
evolution in Q 2 is governed by the DGLAP equations (Gri- 
bov and Lipatov 1972; Lipatov 1975; Altarelli and Parisi, 
1977; Dokshitzer, 1977): 

dDi{x,Q 2 ) = 
d\nQ 2 

^rY, jffjttWCftWf.Q 3 ), CD 

where Pji(z) is the splitting function. We solve the DGLAP 
equations numerically with the conventional QCD (SM case) 
splitting functions and with the supersymmetric (MSSM ca- 
se) ones (Jones and Llewellyn Smith 1983). For the top and 
the MSSM partons we used the FFs of ref. (Rubin 1999). 
While solving the DGLAP equations each parton is included 
at its own threshold energy. 

At small values of x, multiple soft gluon emission can be 
described by the MLLA. It describes the observed hadropro- 
duction quite accurately in the small x region (see eg. Abreu 
et al., 1999; Abbiendi et al., 2000). For large values of x the 
MLLA should not be used. We smoothly connect the solu- 
tion for the FF obtained by the DGLAP equations and the 
MLLA result at a given x c value. Our final result on mx is 
rather insensitive to the choice of x c , the uncertainty is in- 
cluded in our error estimate. We also determined the FF of 
the pion. Fig. [I] shows the FF for the proton and pion at 
Q = 10 16 GeV in SM and MSSM. 
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Figure 1 : The FFs averaged over the quark flavors at Q = 
10 16 GeV for proton/pion in SM (solid/dotted line) and in 
MSSM (dashed/dashed-dotted line) in the relevant x region. 
To show both the small and large x behavior we change from 
logarithmic scale to linear at x = 0.01. 

3 Comparison of the predicted and the 
observed spectra 

UHECR protons produced in the halo of our galaxy can prop- 
agate unaffected and the production spectrum should be com- 
pared with the observations. 

Particles of EG origin and energies above w 5 ■ 10 19 eV 
loose a large fraction of their energies due to interactions 
with CMBR (Greisen, 1966; Zatsepin and Kuzmin, 1966). 
This effect can be quantitatively described by the function 
P(r, E, E c ), the probability that a proton created at a dis- 
tance r with energy E arrives at Earth above the threshold 
energy E c (Bahcall and Waxman 2000). This function has 
been calculated for a wide range of parameters in (Fodor and 
Katz, 2001a), which we use in the present calculation. The 
original UHECR spectrum is changed at least by two differ- 
ent ways: (a) there should be a steepening due to the GZK 
effect; (b) particles loosing their energy are accumulated just 
before the cutoff and produce a bump. We study the ob- 
served spectrum by assuming a uniform source distribution 
for UHECRs. 

Our analysis includes the published and the unpublished 
(from the www pages of the experiments) UHECR data of 
AGASA, Fly's Eye, Haverah Park and HiRes. Since the de- 
cay of SPs results in a non-negligible flux for lower energies 
\og(E m i n / eW) = 18.5 is used as a lower end for the UHECR 
spectrum. Our results are insensitive to the definition of the 
upper end (the flux is extremely small there) for which we 



2 




20 21 22 
Log (E [eV]) 




Figure 2: The available UHECR data with their error bars 
and the best fit from a decaying SP. Note that there are no 
events above 3 x 10 20 eV (shown by an arrow). Nevertheless 
the experiments are sensitive even in this region. Zero event 
does not mean zero flux, but a well defined upper bound for 
the flux (given by the Poisson distribution). Therefore the 
experimental value of the integrated flux is in the "hatched" 
region with 68% confidence level, ("hatching" is a set of 
individual error bars; though most of them are too large to be 
depicted in full) Clearly, the error bars are large enough to be 
consistent with the SP decay. 

choose \og(E max / eV) = 26. As it is usual we divided each 
logarithmic unit into ten bins. Using a Monte-Carlo method 
we included this uncertainty in the final error estimates. The 
predicted number of events in a bin is given by 



N(i) 



/ [A-E 

JEi 



-3.16 



B-j(E,m x )], (2) 



where Ei is the lower bound of the \ th energy bin. The first 
term describes the data below 10 19 eV according to (Takeda 
et al., 1998), where the SP decay gives negligible contribu- 
tion. The second one corresponds to the spectrum of the de- 
caying SPs. A and B are normalization factors. 

The expectation value for the number of events in a bin 
is given by eqn. ([|) and it is Poisson distributed. To deter- 
mine the most probable mx value we used the maximum- 
likelihood method by minimalizing the x 2 (A, B, mx) for 
Poisson distributed data (Groom et al., 2000) 

26.0 

X 2 = 2[N(i)-N (i) + N (i)ln(N (i)/N(i))}, 

i=18.5 

(3) 

where N a (i) is the total number of observed events in the 
\ th bin. In our fitting procedure we have three parameters: 



Figure 3: The most probable values for the mass of the 
decaying ultra heavy dark matter with their error bars and the 
total x 2 values. Note that 21 bins contain nonzero number of 
events and eqn.(Q) has 3 free parameters. 

A, B and mx- Fig. || shows the measured UHECR spectrum 
and the best fit in the EG-MSSM scenario. The first bump 
of the fit represents particles produced at high energies and 
accumulated just above the GZK cutoff due to their energy 
losses. The bump at higher energy is a remnant of mx- In 
the halo models there is no GZK bump, so the relatively large 
x part of the FF moves to the bump around 5 x 10 19 GeV 
resulting in a much smaller mx than in the EG case. The ex- 
perimental data is far more accurately described by the GZK 
effect (dominant feature of the EG fit) than by the FF itself 
(dominant for halo scenarios). 



4 Results 

To determine the most probable value for the mass of the SP 
we studied 4 scenarios. Fig. || contains the Xmm va l ues an d 
the most probable masses with their errors for these scenar- 
ios. 

The UHECR data favors the EG-MSSM scenario. The 
goodnesses of the fits for the halo models are far worse. The 
SM and MSSM cases do not differ significantly. The most 
important message is that the masses of the best fits (EG 
cases) are compatible within the error bars with the MSSM 
gauge coupling unification GUT scale (Amaldi at al., 1991). 

The SP decay will also produce a huge number of pions 
which will decay into photons. Our spectrum contains 94% 
of pions and 6% of protons. This ir/p ratio is in agree- 
ment with (Bhattacharjee and Sigl, 2000; Bhattacharjee and 
Sigl, 2000) which showed that for different classes of mod- 
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els mx < 10 16 GeV , which is the upper boundary of our 
confidence intervals, the generated gamma spectrum is still 
consistent with the observational constraints. We performed 
the whole analysis including the pion produced 7-s in eqn. 
([3]). The results agree with our results of Fig. |3| within error- 
bars. 

In the near future the UHECR statistics will probably be 
increased by an order of magnitude (Boratav, 1996; Guerard, 
1999; Bertou et al. 2000). Performing our analysis for such a 
statistics the uncertainty of mx was found to be reduced by 
two orders of magnitude. 

Since the decay time should be at least the age of the uni- 
verse it might happen that such SPs overdose the universe. 
Due to the large mass of the SPs a single decay results in 
a large number of UHECRs, thus a relatively small number 
of SPs can describe the observations. We checked that in 
all of the four scenarios the minimum density required for 
the best-fit spectrum is more than ten orders of magnitude 
smaller than the critical one. 

More details of the present analysis can be found in (Fodor 
and Katz, 2001b). Note, that a similar study based on the Z- 
burst scenario (Fargion et al., 1999; Weiler, 1999) can be car- 
ried out which gives the mass of the heaviest neutrino (Fodor 
et al., 2001). 
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